Chaotic orbits in a 3D galactic dynamical model with a double nucleus 
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Abstract 



o 

5— 1 . A 3D-dynamical model is constructed for the study of motion in the central regions of a disk galaxy with 
a double nucleus. Using the results of the 2D-model, we find the regions of initial conditions in the 
(x,p x , z,p y ) — Ej, (y = p z = 0) phase space producing regular or chaotic orbits. The majority of stars 
I are on chaotic orbits. All chaotic orbits come arbitrary close to one or to both nuclei. Regular orbits are 
those starting near the stable periodic orbits of the 2D-system with small values of zq. All regular orbits 
i circulate around only one of the two nuclei. 
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O . 1. Introduction 

Astronomers are not sure how common the double nucleus phenomenon is. This happens because the 
centers of galaxies are obscured by dust and gas and, therefore, have yield few of their secrets. Fortunately 
this is beginning to change with the advent of space-based observations. Galaxies with reported double 
nuclei are M31 (Statler et al., 1999), the barred spiral galaxy M83 (Mast et al., 2006; Soria & Wu, 2002), 
NGC 6240 (Komosca et al., 2003; Risaliti et al., 2006), NGC 3256 (Lira et al., 2007) and the dwarf elliptical 
galaxy VCC 128 (Debattista et al., 2006). A detailed catalog of disk galaxies with double nuclei was compiled 
^ ■ by Cimeno et al. (2004). 

On this basis, there is no doubt that the construction of a 3D dynamical model in order to study the 
motion in a galaxy with a double nucleus would be of interest. As far as the authors know there are few 
dynamical models for the study of motion in galaxies with double nuclei (see Jalali & Rafie, 2001; Samlhus 
(N ' & Sridhar, 2002). 

This article can be considered as a continuation of the results presented in a recent paper (see Caranicolas 
& Papadopoulos, 2009; hereafter CP). In Section 2 we present the dynamical model. In Section 3 the regular 
or chaotic character of orbits in the 2D model is explored. In Section 4 we study the nature of motion in 
the 3D model. A discussion and the conclusions of this research are presented in Section 5. 



2. The 3D-model 

The model is an extension of the model used in CP in the 3D space with an additional disk. Thus the 
potential of the first body is 
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where M4, M n \ is the mass of disk and nucleus 1 respectively, a is the disk's scale length, h is the disk's 
scale height, while c„i is the scale length of nucleus 1. The second nucleus is described by the potential 
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where M n2 , c n2 is the mass and the scale length of the nucleus 2 respectively. As in CP, the two bodies 
move in circular orbits in an inertial frame OXYZ with the origin at the center of mass of the system at a 
constant angular velocity ft p > 0. 

Let M t = Md + M n \ + M n2 be the total mass of the system and R be the distance between the two 
bodies. A clockwise rotating frame Oxyz is used with axis Oz coinciding with the axis OZ and the axis Ox 
coinciding with the straight line joining the two bodies. In this frame, which rotates at an angular velocity 
il p , the two bodies have fixed positions £1,2/1 =0 and £2,1/2 = 0. The total potential responsible for the 
motion of a test particle (star) with a unit mass is 
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The angular frequency fl p is calculated as follows: The two bodies move about the center of mass of the 
system with angular velocities Q p2 given by 
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where r 2 = x 2 + y 2 . As the two bodies are not mass points the two angular frequencies are not equal. 
Nevertheless, W6 Ccin niciks th.6m 6qu£l1 by choosing properly the parameters o;, c n i, c n 2 of the system. 
The authors must make clear that, after choosing properly the parameters the two frequencies may differ 
slightly, so that v = (£l p i — fi p2 )/^ P i is of order of 10 -6 or smaller. Therefore, we consider the two angular 
frequencies practically equal, that is fi p i = £l p2 = Q p . 
The equations of motion are written 



d<f> n „ . .. <9$ „^ . .. a$ 
X = ~~dx~ V ^~~dy + pX ' Z = ~~dz' 



(10) 



where the dot indicates derivative with respect to the time. The only integral of motion for the system of 
differential equations (10) is 



J = ~ {Pi +P 2 y + vl) + 2/, z) = Ej, 



(11) 



where p x , p y and p z are the momenta per unit mass conjugate to x, y and z respectively. This is the well 
known Jacobi integral and Ej is its numerical value. In this work we use the same system of galactic units, 
as in CP. The unit of length is 1 kpc, the unit of mass is 2.325 x 10 7 M Q and the unit of time is 0.977 x 10 8 
yr. The velocity unit is equal to 10 km/s, while G is equal to unity. 
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Figure 1: The (x,p x ) phase plane when M d = 1100, M nl = 100, M n2 = 600, R = 1.5, c„i = 0.10, c n2 = 0.35, fl p = 22.1938, a = 
0.3053, h = 0.06. The value of E J2 is equal to -2010. 

3. Orbits in the 2D potential 

In this section we shall study the 2D potential. Thus we set z = p z = in (11) obtaining 

J2 = l{pl+p 2 y )+<S>(x,y) = E J2 , (12) 

where Ej 2 is the numerical value of J 2 . As the system is now two-dimensional, we can use the classical 
method of the (x,p x ), y = 0,p y > 0, Poincare phase plane in order to explore the regular or chaotic character 
of motion. The results obtained from the study of the 2D system will be used in order to determine the 
character of orbits in the 3D system. 

Our results come from the numerical integration of the equations of motion, which was done using a 
Bulirsh-Stoer method in double precision. The accuracy of the calculations was checked by the constancy 
of the Jacobi integral, which was conserved up to the twelfth significant figure. 

Fig. 1 shows the (x,p x ) phase plane when Md = 1100, M n \ — 100, M n2 = 600, R = 1.5, c n \ = 0.10, c n2 = 
0.35, il p — 22.1938, a = 0.3053, h — 0.06, while Ej 2 = —2010. As one can see there is a large unified chaotic 
sea surrounding both nuclei. One also observes two separate regular regions near each nucleus. It is of 
interest to note that the regular regions are around the retrograde periodic points in each of the two nuclei. 
Some smaller regular regions are also present near the smaller nucleus 1 on the left. One additional feature is 
the presence of several small islands near the heavy nucleus 2 on the right. These small islands are produced 
by secondary resonances. 

Fig. 2 is similar to Fig. 1 but when M d = 1100, M nl = 300, M n2 = 400, R = 1.5, c nl = 0.10, c n2 = 
0.25, fip = 22.6244, a = 0.2181, h = 0.06, while Ej 2 = —2100. Here again we see a large unified chaotic sea 
surrounding both nuclei. The two regular areas around the retrograde periodic points of each nucleus are 
also present. On the other hand, a careful observer is able to see some significant differences between the 
two patterns. In Fig. 2 the regular area near the nucleus 1 on the left is now considerably larger than that 
observed in Fig. 1. Furthermore an additional small regular region has appeared near direct periodic point 
of the nucleus 1 on the left. Moreover the corresponding regular area near the nucleus 2 on the right has 
now become smaller. 

This can be explained, if we take into account that the observed chaotic motion near each nucleus is 
a result not only of the force coming from the nucleus itself but also of the force coming from the other 
nucleus. As the mass of nucleus 2 is smaller while the mass of the disk and the distance R between the two 
nuclei are the same, the influence from this nucleus to the disk and the nucleus 1 system is smaller. This 
has as a result the reduction of chaotic region and a parallel increase of the regular regions associated with 
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Figure 2: Similar to Fig. 1 but when M d = 1100, M„i = 300, M n2 = 400, R = 1.5, c„i = 0.10, c n2 = 0.25, U p = 22.6244, a = 
0.2181, h = 0.06, while £ J2 = -2100. 

the nucleus 1 and the disk. On the other hand, for similar reasons, as the mass of nucleus 1 has increased, 
the regular area associated with nucleus 2 has become smaller. 

The numerical results suggest that there are two kinds of chaotic orbits: (i) chaotic orbits approaching 
both nuclei and (ii) chaotic orbits that approach only one of the two nuclei. On the other hand, the regular 
orbits circulate around only one of the two nuclei. 

4. Orbits in the 3D potential 

In this section we shall study the properties of orbits in the 3D potential. In order to keep things simple 
we shall use our experience gained from the study of the 2D system in order to obtain a clear picture of the 
properties of orbits in the 3D dynamical model. 

We are particularly interested to locate the initial conditions in the 3D model producing regular or 
chaotic orbits. A convenient way to obtain this is to start from the (x,p x ) phase plane of the 2D system 
with the same value of the Jacobi integral used in the 2D system. Thus we take Ej = Ej 2 . For this purpose 
a large number of orbits were computed with initial conditions {xo 1 p X Q 1 z a), where {xo,p x o) is a point in the 
chaotic sea of Figs. 1 or 2 with all permissible values of z , and p zn = 0. Remember that, as we are on 
the phase plane, we have yo — 0, while in all cases the value of p y o was obtained from the Jacobi integral. 
All tested orbits were found to be chaotic. Therefore, one concludes that the majority of orbits in the 3D 
system are chaotic. 

An interesting question one might ask is this. Are there any other chaotic orbits in the 3D system? In 
order to give an answer we have taken the sections of the 3D orbit with the plane y = 0, when p y > 0. The 
set of the four-dimensional points (x,p x , z,p z ) was projected on the (z,p z ) plane. If the projected points lie 
on an "invariant curve" this suggests that the motion is regular, while, if not, this is an indication that the 
motion is chaotic. Fig. 3 shows such "invariant curves" for orbits starting near the regular region on the right 
side of Fig. 1. In order to obtain the results shown in Fig. 3 we have taken the point (xo,p x o) = (—0.17, 0.0) 
representing approximately the position of the periodic orbits in the (x,p x ), y = 0,p y > 0, phase plane and a 
set of values of z = (0.02, 0.05, 0.08, 0.11, 0.15, 0.20, 0.25). Note that the numerical results indicate that, for 
small values of z the motion is regular, while for larger values of z , the motion seems to be chaotic. Here 
we must emphasize, that the results of Fig. 3 are rather qualitative and can be considered as an indication 
that the transition from regularity to chaos occurs as the value of zq increases. Results not given here show 
a similar behavior near each regular region in Figs. 1 and 2. 

In order to estimate the degree of chaos in the 2D and the 3D system we have calculated the maximum 
LCE (Lyapunov Characteristic Exponent) (see Lichtenberg & Lieberman, 1992) for a large number of chaotic 
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Figure 3: Projection of the sections of the 3D orbit with the plane y = 0, when p y > 0. The set of the four dimensional points 
(x,p x , z,p z ) is projected on the (z,p z ) plane. 




Figure 4: (a-b): (a-left) A chaotic orbit approaching both nuclei. The initial conditions are xq = 0.8, yo = p x Q = 0,zo = 
0.1, PzO = 0. The values of the other parameters are as in Fig. 1, while Ej = —2010. (b-right) A regular orbit circulating 
around nucleus 1 with initial conditions xo = 1.4, yo = p x Q = 0, zq = 0.15,p z o = 0. The values of the other parameters are as 
in Fig. 2 and Ej = -2100. 

orbits. Each LCE was calculated for a time period of 10 4 time units. The LCE for the 2D orbits, starting 
in the chaotic sea of Fig. 1, was found in the range 5.0—5.2, while in the chaotic sea of Fig. 2 it was found 
in the range 6.0 — 6.2. The LCE for chaotic orbits in the 3D system with initial conditions (x ,p x o, z n ) with 
(xq,Pxo) in the chaotic sea of the Fig. 1 was found in the range 3.2 — 3.4. The corresponding values of LCE 
for Fig. 2 were found in the range 3.9 — 4.1. 

Fig. 4a shows a chaotic orbit with conditions x — 0.8, yo = p x o = 0, z = 0.1, p z o = 0. The value of 
p y o is always found from the Jacobi integral. The values of the other parameters are as in Fig. 1, while 
Ej = —2010. Note that the orbit goes arbitrary close to both nuclei. It is interesting to observe that near 
the more massive nucleus the orbit is deflected to more higher values of z, while near the less massive nucleus 
the star stays close to the disk. Fig. 4b shows a regular orbit circulating around nucleus 1. The initial 
conditions are x = 1.4, y = Pxo = 0, z = 0.15,p zu = 0. The values of the other parameters are as in Fig. 



Figure 5: (a-b): (a-left) A regular orbit and (b-right) a chaotic orbit. The two orbits differ in initial conditions only in the 
value of zq. See text for details. 




2, while Ej = -2100. 

Fig. 5a shows a quasi-periodic orbit starting near the retrograde periodic point, which is close to nucleus 
1. The initial conditions are x = 0, yo = Pxo = 0, z = 0.1, p z o = 0. The values of the other parameters are 
as in Fig. 2, while Ej = —2100. Fig. 5b shows an orbit with the same initial conditions, the same value 
of the Jacobi integral and the same values of the parameters as in Fig. 5a but when zq = 0.5. The orbit 
has now become chaotic and goes arbitrarily close to the nucleus 1. This orbit shows, from another point 
of view, that 3D orbits starting near the stable periodic points of the 2D system are regular only for small 
values of zq ■ 

Here the physical parameter playing an important role is the L z component of the angular momentum. 
From our previous experience we know that low angular momentum stars, on approaching a dense nucleus 
are scattered off the galactic plane displaying chaotic motion (Caranicolas & Innanen, 1991; Caranicolas & 
Papadopoulos, 2003). Of course here in 3D space things are more complicated than in an axially symmetric 
dynamical model, where L z is conserved (see Caranicolas & Innanen, 1991). As L z , is not conserved, we 
can compute numerically the mean value < L z > of L z using the formula 

1 " 

<L z >=-Y j L zi . (13) 

i 

Our numerical calculations suggest that the chaotic orbits have low values of < L z >, while regular 
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orbits have high values of < L z >. Fig. 6a shows the evolution of L z with the time for the regular orbit of 
Fig. 5a. Here we find a value of < L z >= 12.5. Fig. 6b is similar to Fig. 6a but for the chaotic orbit of 
Fig. 5b. Here < L z >= —10.3 The value of n in both cases was 10 4 . 

5. Discussion and conclusions 

Observation data show that a small fraction of active galaxies have double nuclei (Eracleus & Halpen, 
2003; Xinwu & Ting-Gui, 2006). It was this reason that motivated our construction of a 3D model in order 
to study the motion in a disk galaxy hosting a binary nucleus. It was found that the majority of orbits in the 
2D system were chaotic. Two kinds of chaotic orbits were observed: (i) Chaotic orbits that approach both 
nuclei and (ii) Chaotic orbits that approach only one of the nuclei. The regular regions are confined mainly 
around the retrograde periodic points in both nuclei. All regular orbits go around nucleus 1 or nucleus 2 
but not both. It was also found that the total velocity near each nucleus attains high values. The value of 
velocity depends on the mass of the nucleus and the value of its scale length. Regular motion corresponds 
to low velocities while chaotic motion is characterized by high velocities. 

In order to understand the behavior of orbits in the 3D system we have used our knowledge from 
the 2D system. Of particular interest was the determination of the region of initial conditions in the 
(x, p Xj z,p y ) = Ej, (y = p z = 0) phase space that produces regular or chaotic 3D orbits. As p y o was found 
always from the Jacobi integral we have used the same value of Ej, as in the 2D system and took initial 
condition (xq,p x q, zq) such as (xq,p x o) lies in the chaotic region of the 2D system. It was found that the 
motion was chaotic for all permissible values of zq. In the case when (xo,p x o) was inside a regular region, 
the corresponding 3D orbit was regular only for small values of z 0} while for larger values of z the orbit 
was chaotic. The particular values of zq were different for each regular region of the 2D system. We did not 
feel that it was necessary to try to define the particular values of zq for each case. 

An important role is played by the L z component of the test particle's angular momentum. It was found 
that the values of < L z > for regular orbits are larger, than those corresponding to chaotic orbits. Thus, the 
L z component of the angular momentum is a significant parameter connected with the regular or chaotic 
character of orbits in both 2D and 3D galactic models. 

In order to estimate the degree of chaos in the 2D as well as in the 3D potentials, we have computed the 
maximum Lyapunov Characteristic Exponent (LCE) for a large number of orbits for a time period of 10 4 
time units. The numerical results indicate that the degree of chaos in 3D double nucleus systems is smaller 
than in similar 2D systems. 
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